#### Load necessary packages ####

install.packages("pacman")

pacman::p_load(tidyverse, magrittr, lubridate, plm, fixest, data.table, readxl, countrycode, texreg, xtable, ggsci, factoextra, FactoMineR, ggcorrplot, corrr,
               cowplot, googleway, ggrepel, ggspatial, sf, lwgeom, rnaturalearth, rnaturalearthdata, kableExtra, Rilostat, httr, gmm, klaR)

#### Set working directory ####

setwd('~') 

#### Load WorkR data for Figures and main data for OLS results ####

workr <- read.csv("WorkR.csv", encoding = "latin1")

OlsData <- read.csv("OlsData.csv")

GmmData <- read.csv("GmmData.csv")

#### Figure 1 ####

workr %>%
  mutate(eu = ifelse((grepl("Austria|Belgium|Bulgaria|Croatia|Republic of Cyprus|Czech Republic|Denmark|Estonia|Finland|France|Germany|Greece|Hungary|Ireland|Italy|Latvia|Lithuania|Malta|Norway|Iceland|Lichtenstein|Poland|Portugal|Romania|Slovakia|Slovenia|Spain|Sweden|Switzerland", country, ignore.case = T))==T, "1", "0"),
         country = ifelse(eu =="1", "Social Europe", country),
         country = ifelse(cowcode %in% OlsData$ccode, "Developing Countries", country)) %>%
  filter(country %in% c("Social Europe", "United States", "Developing Countries")) %>%
  dplyr::select(country, year, uri_p, sri_p, uri_l, sri_l) %>%
  mutate(uri = (uri_p+uri_l)/2,
         sri = (sri_p+sri_l)/2) %>%
  pivot_longer(uri:sri, names_to = "Right", values_to = "Value") %>%
  mutate(Right = ifelse(Right == "uri", "Union Rights", "Substantive Rights")) %>%
  ggplot() + 
  stat_summary(aes(x = year, y = Value, group = country, color = country), fun = "mean", geom = "smooth") +
  theme_minimal() +
  labs(x = "Year", y = "WorkR Score", color = "") +
  theme(axis.text.x = element_text(angle = 45)) +
  facet_wrap(~factor(Right, levels = c('Union Rights', 'Substantive Rights')), scales = "free")

#### Figure 2 ####
### Produces some warnings regarding unmatched country codes and ###
### country names. None of these are included in the analysis.  ###

d <- workr %>%
  group_by(country) %>%
  summarize(uri_p = mean(uri_p, na.rm = T),
            sri_p = mean(sri_p, na.rm = T)) %>%
  ungroup() %>%
  mutate(ccode = as.character(countrycode(country, "country.name", "cown")),
         region = as.character(countrycode(as.numeric(ccode), "cown", "un.regionsub.name"))) %>%
  dplyr::select(-country) %>%
  filter(is.na(ccode)==F,
         ccode %in% OlsData$ccode)

world <- ne_countries(scale = "medium", returnclass = "sf")

world2 <- world %>%
  mutate(ccode = as.character(countrycode(name_long, 'country.name', 'cown'))) %>%
  left_join(d, by = 'ccode') 

ggpubr::ggarrange(
  ggplot(data = world2) +
    geom_sf(aes(fill = uri_p), size = 0.25) + 
    theme_bw() +
    scale_fill_gradient(low = 'white', high= '#194771') +
    theme(panel.background = element_rect(fill = 'white'),
          legend.position = 'right') + 
    labs(fill = "Union Rights"),
  ggplot(data = world2) +
    geom_sf(aes(fill = sri_p), size = 0.25) + 
    theme_bw() +
    scale_fill_gradient(low = 'white', high= '#194771') +
    theme(panel.background = element_rect(fill = 'white'),
          legend.position = 'right') + 
    labs(fill = "Substantive Rights"),
  nrow = 2,
  align = "v"
)

#### Figure 3 ####

OlsData %>%
  mutate(region2 = case_when(grepl("Eastern Asia|South-eastern Asia|Polynesia|Melanesia", region) ~ "East Asia",
                             grepl("Western Asia|Central Asia|Southern Asia", region) ~ "South Asia",
                             grepl("Europe", region) ~ "Eastern Europe",
                             TRUE ~ region)) %>%
  group_by(region2, year) %>%
  summarize(across(c(uri_p,sri_p), \(x) mean(x, na.rm = T))) %>%
  ungroup() %>%
  filter(complete.cases(.)) %>%
  pivot_longer(uri_p:sri_p, names_to = "Right", values_to = "Value") %>%
  mutate(Right = ifelse(Right == "uri_p", "Union Rights", "Substantive Rights")) %>%
  ggplot() +
  geom_smooth(aes(x = year, y = Value, group = region2, color = region2, fill = region2), linewidth = 0.75, se = F) +
  theme_minimal() +
  labs(color = "", fill = "") +
  facet_wrap(~factor(Right, levels = c('Union Rights', 'Substantive Rights')), scales = "free")

#### Figure 4 ####

OlsData %>%
  dplyr::select(year, eu_stock, usa_stock) %>%
  group_by(year) %>%
  mutate(across(contains("stock"), \(x) sum(x, na.rm = T))) %>%
  ungroup() %>%
  distinct() %>%
  pivot_longer(-year, names_to = "Host", values_to = "Stock FDI") %>%
  ggplot() +
  geom_smooth(aes(x = year, y = `Stock FDI`, group = Host, color = Host), linewidth = 0.75, method = "loess", se = F) +
  theme_minimal() +
  scale_color_discrete(name = "Host", labels = c("Social Europe", "USA")) + 
  labs(x = "Year", y = "FDI Stock (Billions)")

#### Figure 5 ####
### Produces some warnings relating to two unmatched country                  ###
### codes and rows automatically removed by ggplot2.                         ###

ggpubr::ggarrange(
  GmmData %>% 
    mutate(region = countrycode(as.numeric(ccode), 'cown', 'un.regionsub.name')) %>%
    filter(is.na(region)==F & region!='Polynesia' & region!='Melanesia') %>%
    mutate(region = ifelse(region == 'Latin America and the Caribbean', 'Latin America and \n the Caribbean', region)) %>%
    ggplot() +
    stat_summary(aes(x = as.character(year), y = GrowthDiff, group = region, color = region),fun = 'mean', geom = 'point', se = F) +
    stat_summary(aes(x = as.character(year), y = GrowthDiff, group = 1), fun = 'mean', geom = 'smooth', color = 'black') +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = 'none') +
    scale_color_jco() +
    xlab('') +
    ylab('Growth Rate Difference (Home Country - Europe)'),
  GmmData %>% 
    mutate(region = countrycode(as.numeric(ccode), 'cown', 'un.regionsub.name')) %>%
    filter(is.na(region)==F & region!='Polynesia' & region!='Melanesia',
           ccode %in% unique(OlsData$ccode)) %>%
    mutate(region = ifelse(region == 'Latin America and the Caribbean', 'Latin America and \n the Caribbean', region)) %>%
    ggplot(aes(x = region, y = ShareDevIfdi, color = region)) +
    stat_summary(fun = 'mean', geom = 'point', group = 1) +
    stat_summary(fun.data = 'mean_cl_normal', geom = 'errorbar', group = 1, width = 0.05, fun.args = list(mult = 1)) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    theme_bw() +
    coord_cartesian(ylim = c(0, 1.0)) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = 'none') +
    scale_color_jco() +
    xlab('') +
    ylab('Share of Developed IFDI (Regional Averages)'),
  GmmData %>%
    mutate(region = countrycode(as.numeric(ccode), 'cown', 'un.regionsub.name')) %>%
    filter(is.na(region)==F & region!='Polynesia' & region!='Melanesia') %>%
    mutate(region = ifelse(region == 'Latin America and the Caribbean', 'Latin America and \n the Caribbean', region)) %>%
    ggplot(aes(x = as.character(year), y = uri_p)) +
    stat_summary(aes(group = region, color = region), fun = 'mean', geom = 'point') +
    stat_summary(aes(group = 1), fun = 'mean', geom = 'smooth', color = 'black') +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = 'none') +
    scale_color_jco() +
    xlab('') +
    ylab('Union Rights in Practice (Regional Averages)'),
  GmmData %>%
    mutate(region = countrycode(as.numeric(ccode), 'cown', 'un.regionsub.name')) %>%
    filter(is.na(region)==F & region!='Polynesia' & region!='Melanesia') %>%
    mutate(region = ifelse(region == 'Latin America and the Caribbean', 'Latin America and \n the Caribbean', region)) %>%
    ggplot(aes(x = as.character(year), y = sri_p)) +
    stat_summary(aes(group = region, color = region), fun = 'mean', geom = 'point') +
    stat_summary(aes(group = 1), fun = 'mean', geom = 'smooth', color = 'black') +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = 'none') +
    scale_color_jco() +
    xlab('') +
    ylab('Substantive Rights in Practice (Regional Averages)')
)

#### Figure 6 ####
### Estimate models first; then the coefficients and other statistics are    ###
### loaded into a data frame which is used to generate the coefficient plot. ###

ivs_vec <- c('eu_stock',
             'eu_stock_lag1',
             'eu_stock_lag2',
             'usa_stock',
             'usa_stock_lag1',
             'usa_stock_lag2')
dvs_vec <- c('assoc_p', 
             'barg_p',
             'force_p',
             'age_p',
             'wage_p',
             'safe_p',
             'hour_p')

ivs <- paste0(" ~ ", ivs_vec)

dvs_ivs <- c(unlist(lapply(ivs, function(x) paste0(dvs_vec, x,
                                                   '+
               v2x_frassoc_thick_lag1 +
               gsp_lag1 + 
               gdp_growth_lag1 +
               UrbanPop_lag1+
               log(1 + ilo_rat_lag1) +
               log(1 + gdppc_wanlin_lag1)+
               log(1 + export_devd_gdp_lag1) + 
               log(1 + InwardFDIstock_lag1) +
               conflict|ccode+year'))))
formulas <- lapply(dvs_ivs, formula)

FirstStupidModels <- lapply(formulas, function(x) {
  feols(x, data = OlsData,
        se = 'cluster')
})


PlotData <- data.frame(outcome = NA, fdi = NA, coefs = NA, se = NA, pvalue = NA, min = NA, max = NA)

PlotSeq <- FirstStupidModels

for(i in seq_along(PlotSeq)) {
  PlotData[i,1] <- as.character(PlotSeq[[i]]$fml[2])
  PlotData[i,2] <- as.character(PlotSeq[[i]][['fml']][[3]][[2]][[2]][[2]][[2]][[2]][[2]][[2]][[2]][[2]])
  PlotData[i,3] <- as.numeric((coef(PlotSeq[[i]])[1]))
  PlotData[i,4] <- as.numeric(PlotSeq[[i]]$se[[1]])
  PlotData[i,5] <- PlotSeq[[i]]$coeftable[1,4]
  PlotData[i,6] <- PlotData[i,3] - PlotData[i,4]
  PlotData[i,7] <- PlotData[i,3] + PlotData[i,4]
}


level_order <- c('t', 't-1', 't-2')

ggpubr::ggarrange(
  PlotData %>%
    filter(outcome %in% c('assoc_p', 
                          'barg_p',
                          'force_p',
                          'age_p',
                          'wage_p',
                          'safe_p',
                          'hour_p')) %>%
    mutate(outcome = case_when(outcome == "assoc_p" ~ "Associational",
                               outcome == "barg_p" ~ "Bargaining",
                               outcome == "force_p" ~ "Forced Labor",
                               outcome == "age_p" ~ "Child Labor",
                               outcome == "wage_p" ~ "Wage Regulations",
                               outcome == "safe_p" ~ "Work Safety", 
                               outcome == "hour_p" ~ "Working Time"),
           fdi2 = str_replace(fdi, 'eu_|usa_', ''),
           Lag = ifelse(fdi2 == 'stock', 't', str_c('t-', str_sub(fdi2, -1)))) %>%
    filter(fdi %in% c('eu_stock', 'eu_stock_lag1', 'eu_stock_lag2', 'eu_stock_lag3')) %>%
    mutate(max95 = coefs+1.96*se,
           min95 = coefs-1.96*se) %>%
    filter(is.na(outcome)==F) %>%
    ggplot() +
    geom_point(aes(y = factor(outcome, 
                              level = c('Forced Labor',
                                        'Child Labor',
                                        'Wage Regulations',
                                        'Work Safety',
                                        'Working Time',
                                        'Bargaining',
                                        'Associational')), x = coefs, group = Lag, color = Lag),
               position = position_dodge(width = 0.2)) +
    geom_errorbar(aes(y = factor(outcome, 
                                 level = c('Forced Labor',
                                           'Child Labor',
                                           'Wage Regulations',
                                           'Work Safety',
                                           'Working Time',
                                           'Bargaining',
                                           'Associational')), xmin = min95, xmax = max95, group = Lag, color = Lag), 
                  width = 0.03, size = 0.2, alpha = 0.75,
                  position = position_dodge(width = 0.2)) +
    geom_vline(xintercept = 0, color = 'black', linetype = 'dashed', alpha = 0.5) +
    coord_cartesian(
      xlim = c(-0.06, 0.06)
    ) +
    theme_classic() +
    theme(legend.position = 'right',
          plot.title = element_text(hjust = 0.5, size = 10)) +
    scale_color_nejm(name = 'Lags') +
    # scale_color_manual(values = c('blue', 'dark blue', 'dark red','dark green'), 
    #                    name = 'Lags') +
    xlab('Coefficients') +
    ylab('Labor Rights') +
    ggtitle('OFDI to EU'),
  PlotData %>%
    filter(outcome %in% c('assoc_p', 
                          'barg_p',
                          'force_p',
                          'age_p',
                          'wage_p',
                          'safe_p',
                          'hour_p')) %>%
    mutate(outcome = case_when(outcome == "assoc_p" ~ "Associational",
                               outcome == "barg_p" ~ "Bargaining",
                               outcome == "force_p" ~ "Forced Labor",
                               outcome == "age_p" ~ "Child Labor",
                               outcome == "wage_p" ~ "Wage Regulations",
                               outcome == "safe_p" ~ "Work Safety", 
                               outcome == "hour_p" ~ "Working Time"),
           fdi2 = str_replace(fdi, 'eu_|usa_', ''),
           Lag = ifelse(fdi2 == 'stock', 't', str_c('t-', str_sub(fdi2, -1))),
           Lag = ifelse(Lag == 't-d', 't+1', Lag)) %>%
    filter(fdi %in% c('usa_stock', 'usa_stock_lag1', 'usa_stock_lag2', 'usa_stock_lag3'),
           is.na(outcome)==F) %>%
    mutate(max95 = coefs+1.96*se,
           min95 = coefs-1.96*se) %>%
    ggplot() +
    geom_point(aes(y = factor(outcome, 
                              level = c('Forced Labor',
                                        'Child Labor',
                                        'Wage Regulations',
                                        'Work Safety',
                                        'Working Time',
                                        'Bargaining',
                                        'Associational')), x = coefs, group = Lag, color = Lag),
               position = position_dodge(width = 0.2)) +
    geom_errorbar(aes(y = factor(outcome, 
                                 level = c('Forced Labor',
                                           'Child Labor',
                                           'Wage Regulations',
                                           'Work Safety',
                                           'Working Time',
                                           'Bargaining',
                                           'Associational')), xmin = min95, xmax = max95, group = Lag, color = Lag), 
                  width = 0.03, size = 0.2, alpha = 0.75,
                  position = position_dodge(width = 0.2)) +
    geom_vline(xintercept = 0, color = 'black', linetype = 'dashed', alpha = 0.5) +
    coord_cartesian(
      # xlim = c(-0.06, 0.06)
    ) +
    theme_classic() +
    theme(legend.position = 'right',
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(),
          axis.title.y =element_blank(),
          plot.title = element_text(hjust = 0.5, size = 10),
          axis.line.y =element_blank()) +
    scale_color_nejm(name = 'Lags') +
    # scale_color_manual(values = c('blue', 'dark blue', 'dark red','dark green'), 
    #                    name = 'Lags') +
    xlab('Coefficients') +
    ylab('Labor Rights') +
    ggtitle('OFDI to USA'),
  common.legend = T,
  legend = 'right',
  widths = c(1.15,1)
)

#### Table 1 ####

ivs_vec <- c('eu_stock_lag1',
             'usa_stock_lag1')
dvs_vec <- c('uri_p',
             'sri_p')

ivs <- paste0(" ~ ", ivs_vec)

dvs_ivs <- c(unlist(lapply(ivs, function(x) paste0(dvs_vec, x,
                                                   '|ccode+year'))),
             unlist(lapply(ivs, function(x) paste0(dvs_vec, x,
                                                   '+
               v2x_frassoc_thick_lag1 +
               gsp_lag1 + 
               log(1 + ilo_rat_lag1) +
               log(1 + gdppc_wanlin_lag1)+
               gdp_growth_lag1 +
               UrbanPop_lag1+
               log(1 + export_devd_gdp_lag1) + 
               conflict +
               log(1 + TotalOfdi_lag1) +
               log(1 + io_mem_n) + 
               log(1 + ngoc)|ccode+year'))),
             unlist(lapply(ivs, function(x) paste0(dvs_vec, x,
                                                   '+
               tfp_lag1 +
               plp_lag1 +
               indva_lag1 +
              gov1rlc_lag1 +
                log(1 + InwardFDIstock_lag1) +
                v2x_frassoc_thick_lag1 +
               gsp_lag1 + 
               log(1 + ilo_rat_lag1) +
               log(1 + gdppc_wanlin_lag1)+
               gdp_growth_lag1 +
               UrbanPop_lag1+
               log(1 + export_devd_gdp_lag1)+
                conflict +
                log(1 + TotalOfdi_lag1) +
                log(1 + io_mem_n) + 
                log(1 + ngoc)|ccode+year'))))
formulas <- lapply(dvs_ivs, formula)

Tables1andA3 <- lapply(formulas, function(x) {
  feols(x, data = OlsData,
        se = 'cluster')
})

screenreg(Tables1andA3[c(1:8)],
        digits = 3,
        custom.header = list('FDI to EU' = 1:2, 'FDI to USA' = 3:4,'FDI to EU' = 5:6, 'FDI to USA' = 7:8),
        custom.model.names = rep(c('URI',
                                   'SRI'),4),
        custom.coef.map = list('eu_stock_lag1'='OFDI Flows t-1',
                               'usa_stock_lag1'='OFDI Flows t-1',
                               'v2x_frassoc_thick_lag1'='Freedom of Assoc. t-1',
                               'gsp_lag1'='GSP t-1',
                               'log(1 + ilo_rat_lag1)'='ILO Ratification t-1',
                               'log(1 + gdppc_wanlin_lag1)'='GDP per capita t-1',
                               'gdp_growth_lag1'= 'GDP Growth t-1',
                               'UrbanPop_lag1'='Urban Population t-1',
                               'log(1 + export_devd_gdp_lag1)'='Exports to Devd t-1',
                               'log(1 + TotalOfdi_lag1)' = 'Total OFDI Stock',
                               'conflict1' = 'Conflict - PRIO',
                               'log(1 + io_mem_n)' = 'IO Membership',
                               'log(1 + ngoc)' = 'INGO Presence'),
        custom.gof.rows = list('Country Fixed Effects' = c('Yes', 'Yes', 'Yes', 'Yes','Yes', 'Yes', 'Yes', 'Yes'),
                               'Year Fixed Effects' = c('Yes', 'Yes', 'Yes', 'Yes','Yes', 'Yes', 'Yes', 'Yes')),
        stars = c(0.01, 0.05, 0.1))

#### Table 2 ####

ivs_vec <- c('ofdi_stock_context_index_lag1',
             'ofdi_stock_context_uri_p_lag1',
             'ofdi_stock_context_sri_p_lag1')
dvs_vec <- c('uri_p',
             'sri_p')

ivs <- paste0(" ~ ", ivs_vec)

dvs_ivs <- c(unlist(lapply(ivs, function(x) paste0(dvs_vec, x,
                                                   '+
               v2x_frassoc_thick_lag1 +
               gsp_lag1 + 
               log(1 + ilo_rat_lag1) +
               log(1 + gdppc_wanlin_lag1)+
               gdp_growth_lag1 +
               UrbanPop_lag1+
               log(1 + export_devd_gdp_lag1) + 
               conflict +
               log(1 + TotalOfdi_lag1) +
               log(1 + io_mem_n) + 
               log(1 + ngoc)|ccode+year'))))
formulas <- lapply(dvs_ivs, formula)

Table2 <- lapply(formulas, function(x) {
  feols(x, data = OlsData,
        se = 'cluster')
})

screenreg(Table2[c(3,6,1,2)],
        digits = 3,
        custom.header = list('De facto labor rights (WorkR dataset)' = 1:2, 'Institution Index (industrial relations)' = 3:4),
        custom.model.names = c('URI',
                               'SRI',
                               'URI',
                               'SRI'),
        custom.coef.map = list('ofdi_stock_context_index_lag1'='OFDI Context t-1',
                               'ofdi_stock_context_uri_p_lag1'='OFDI Context t-1',
                               'ofdi_stock_context_sri_p_lag1'='OFDI Context t-1',
                               'v2x_frassoc_thick_lag1'='Freedom of Assoc. t-1',
                               'gsp_lag1'='GSP t-1',
                               'log(1 + ilo_rat_lag1)'='ILO Ratification t-1',
                               'log(1 + gdppc_wanlin_lag1)'='GDP per capita t-1',
                               'gdp_growth_lag1'= 'GDP Growth t-1',
                               'UrbanPop_lag1'='Urban Population t-1',
                               'log(1 + export_devd_gdp_lag1)'='Exports to Devd t-1',
                               'log(1 + TotalOfdi_lag1)' = 'Total OFDI Stock',
                               'conflict1' = 'Conflict - PRIO',
                               'log(1 + io_mem_n)' = 'IO Membership',
                               'log(1 + ngoc)' = 'INGO Presence'),
        custom.gof.rows = list('Country Fixed Effects' = c('Yes', 'Yes','Yes', 'Yes'),
                               'Year Fixed Effects' = c('Yes', 'Yes','Yes', 'Yes')),
        stars = c(0.01, 0.05, 0.1)
)

#### Table 3 ####
### The estimates for GMM models 3 and 4 produce warnings regarding their    ###
### default solution to the singularity of the first-step matrix, which is   ### 
### to use a general inverse.                                                ###    

pgmm1 <- plm::pgmm(uri_p ~ plm::lag(eu_stock, 1) +
                     plm::lag(uri_p, 1:2) +
                     plm::lag(v2x_frassoc_thick,1)+
                     plm::lag(UrbanPop, 1) +
                     plm::lag(gdp_growth, 1) +
                     plm::lag(gsp, 1) +
                     plm::lag(log(1 + gdppc_wanlin),1) +
                     plm::lag(log(1 + export_devd_gdp), 1) +
                     plm::lag(log(1 + ilo_rat), 1) +
                     plm::lag(log(1 + ngoc), 1) +
                     plm::lag(as.numeric(conflict), 1) +
                     plm::lag(log(1 + io_mem_n), 1) +
                     plm::lag(log(1 + TotalOfdi), 1)
                   |plm::lag(uri_p, 3:99) + plm::lag(eu_stock,2:99)
                   |plm::lag(I(GrowthDiff*ShareDevIfdi), 2:8)
                   , data = GmmData, 
                   collapse = T,
                   model = 'twostep',
                   effect = 'twoway', 
                   index = c('ccode', 'year'),
                   transform = 'ld'
)

s1 <- summary(pgmm1, robust = T)

pgmm2 <- plm::pgmm(sri_p ~ plm::lag(eu_stock, 1) +
                     plm::lag(sri_p, 1:2) +
                     plm::lag(v2x_frassoc_thick,1)+
                     plm::lag(UrbanPop, 1) +
                     plm::lag(gdp_growth, 1) +
                     plm::lag(gsp, 1) +
                     plm::lag(log(1 + gdppc_wanlin),1) +
                     plm::lag(log(1 + export_devd_gdp), 1) +
                     plm::lag(log(1 + ilo_rat), 1) +
                     plm::lag(log(1 + ngoc), 1) +
                     plm::lag(as.numeric(conflict), 1) +
                     plm::lag(log(1 + io_mem_n), 1) +
                     plm::lag(log(1 + TotalOfdi), 1)
                   |plm::lag(sri_p, 3:99) + plm::lag(eu_stock,2:99)
                   |plm::lag(I(GrowthDiff*ShareDevIfdi), 2:8)
                   , data = GmmData, 
                   collapse = T,
                   model = 'twostep',
                   effect = 'twoways', 
                   index = c('ccode', 'year'),
                   transform = 'ld'
)


s2 <- summary(pgmm2, robust = T)

pgmm3 <- plm::pgmm(uri_p ~ plm::lag(usa_stock, 1) +
                     plm::lag(uri_p, 1:2) +
                     plm::lag(v2x_frassoc_thick,1)+
                     plm::lag(log(1 + UrbanPop), 1) +
                     plm::lag(gdp_growth, 1) +
                     plm::lag(gsp, 1) +
                     plm::lag(log(1 + gdppc_wanlin),1) +
                     plm::lag(log(1 + export_devd_gdp), 1) +
                     plm::lag(log(1 + ilo_rat), 1) +
                     plm::lag(log(1 + ngoc), 1) +
                     plm::lag(as.numeric(conflict), 1) +
                     plm::lag(log(1 + io_mem_n), 1) +
                     plm::lag(log(1 + TotalOfdi), 1)
                   |plm::lag(uri_p, 3:99) + plm::lag(usa_stock,2:99)
                   |plm::lag(I(GrowthDiff*ShareDevIfdi), 2:8)
                   , data = GmmData, 
                   collapse = T,
                   model = 'twostep',
                   effect = 'twoways', 
                   index = c('ccode', 'year'),
                   transform = 'ld'
)

s3 <- summary(pgmm3, robust = T)

pgmm4 <- plm::pgmm(sri_p ~ plm::lag(usa_stock, 1) +
                     plm::lag(sri_p, 1:2) +
                     plm::lag(v2x_frassoc_thick,1)+
                     plm::lag(log(1 + UrbanPop), 1) +
                     plm::lag(gdp_growth, 1) +
                     plm::lag(gsp, 1) +
                     plm::lag(log(1 + gdppc_wanlin),1) +
                     plm::lag(log(1 + export_devd_gdp), 1) +
                     plm::lag(log(1 + ilo_rat), 1) +
                     plm::lag(log(1 + ngoc), 1) +
                     plm::lag(as.numeric(conflict), 1) +
                     plm::lag(log(1 + io_mem_n), 1) +
                     plm::lag(log(1 + TotalOfdi), 1)
                   |plm::lag(sri_p, 3:99) + plm::lag(usa_stock,2:99)
                   |plm::lag(I(GrowthDiff*ShareDevIfdi), 2:8)
                   , data = GmmData, 
                   collapse = T,
                   model = 'twostep',
                   effect = 'twoways', 
                   index = c('ccode', 'year'),
                   transform = 'ld'
)

s4 <- summary(pgmm4, robust = T)

screenreg(list(pgmm1,pgmm2,pgmm3,pgmm4),
        override.se = list(s1$coefficients[,2],
                           s2$coefficients[,2],
                           s3$coefficients[,2],
                           s4$coefficients[,2]),
        override.pvalues = list(s1$coefficients[,4],
                                s2$coefficients[,4],
                                s3$coefficients[,4],
                                s4$coefficients[,4]),
        custom.header = list('FDI to EU' = 1:2, 'FDI to USA' = 3:4),
        custom.model.names = c('Union Rights', 'Substantive Rights', 'Unions Rights', 'Substantive Rights'),
        custom.coef.map = list('plm::lag(eu_stock, 1)'='OFDI Stockt-1',
                               'plm::lag(usa_stock, 1)'='OFDI Stockt-1'),
        stars = c(0.01, 0.05, 0.1),
        digits = 3,
        custom.gof.rows = list('All Controls' = c('Yes', 'Yes', 'Yes', 'Yes'),
                               'Country Fixed Effects' = c('Yes', 'Yes', 'Yes', 'Yes'),
                               'Year Fixed Effects' = c('Yes', 'Yes', 'Yes', 'Yes'),
                               'Arellano-Bond test for AR(1)'= c(s1$m1$statistic,
                                                                 s2$m1$statistic,
                                                                 s3$m1$statistic,
                                                                 s4$m1$statistic),
                               'Arellano-Bond test for AR(2)'= c(s1$m2$statistic,
                                                                 s2$m2$statistic,
                                                                 s3$m2$statistic,
                                                                 s4$m2$statistic))
)

#### Figure A1 ####
### First estimate the base model and extract the residuals which are used ###
### as weights. The weights are then joined back to the main data which is ###
### then joined with the loaded world map. There are some locations        ###
### (primarily micro-states and territories) that do not match with the    ###
### correlates of war codes (cown).                                        ###  

RegWt <- feols(uri_p ~
                 v2x_frassoc_thick_lag1 +
                 gsp_lag1 + 
                 gdp_growth_lag1 +
                 UrbanPop_lag1+
                 log(1 + ilo_rat_lag1) +
                 log(1 + gdppc_wanlin_lag1)+
                 log(1 + export_devd_gdp_lag1) + 
                 log(1 + InwardFDIstock_lag1) +
                 conflict|ccode+year,
               data = OlsData,
               se = 'hetero')

Wts <- RegWt$residuals
Obsvs <- c(1:1464)
ObsvsRm <- abs(c(RegWt$obs_selection$obsRemoved))
Obsvs <- Obsvs[-ObsvsRm]

Wts <- data.frame(w = Wts, Obs = as.character(Obsvs))

OlsData$Obs <- row.names(OlsData)

OlsData %<>%
  left_join(Wts)

CountryWs <- OlsData %>%
  mutate(CountryName = countrycode(as.numeric(ccode), 'cown', 'country.name'),
         w = w^2) %>%
  group_by(ccode) %>%
  mutate(CountryWa = sum(w, na.rm = T),
         DeltaEuStock = max(eu_stock)-min(eu_stock)) %>%
  ungroup() %>%
  dplyr::select(ccode, CountryName, year, w, CountryWa, eu_stock, DeltaEuStock) %>%
  mutate(AllW = sum(w, na.rm = T),
         CountryW = CountryWa/AllW) %>%
  distinct(CountryName, .keep_all = T)

world <- ne_countries(scale = "medium", returnclass = "sf")

world2 <- world %>%
  mutate(ccode = countrycode(iso_a3, 'iso3c', 'cown')) %>%
  left_join(CountryWs, by = 'ccode') %>%
  mutate(CtwBi = as.character(ifelse(CountryW > 0, 1, 0)),
         CountryWInt = CountryW*100)

ggplot(data = world2) +
  geom_sf(aes(fill = CountryWInt), size = 0.25) + 
  theme_bw() +
  scale_fill_gradient(low = 'white', high= '#194771', limits = c(0, 3)) +
  theme(panel.background = element_rect(fill = 'white'),
        legend.position = 'none')

#### Figure A2 ####
### First, run regressions iteratively omitting each country. Coefficients ###
### are then extracted and plotted in the histogram.                       ###

countries <- unique(OlsData$ccode)

CtrySubsetResults <- data.frame(OmittedCountry = rep(NA, 122), Coef = rep(NA, 122), SE = rep(NA, 122), TStat = rep(NA, 122), PVal = rep(NA, 122))

for(i in seq_along(countries)) {
  fereg <- feols(uri_p ~
                   eu_stock_lag1 +
                   v2x_frassoc_thick_lag1 +
                   gsp_lag1 + 
                   gdp_growth_lag1 +
                   UrbanPop_lag1+
                   log(1 + ilo_rat_lag1) +
                   log(1 + gdppc_wanlin_lag1)+
                   log(1 + export_devd_gdp_lag1) + 
                   log(1 + InwardFDIstock_lag1) +
                   conflict|ccode+year,
                 data = OlsData,
                 subset = (OlsData$ccode!=countries[i]),
                 se = 'cluster')
  CtrySubsetResults[i,1] <- countrycode(as.numeric(countries)[i], 'cown', 'country.name')
  CtrySubsetResults[i,2] <- as.numeric(fereg$coeftable[1,1])
  CtrySubsetResults[i,3] <- as.numeric(fereg$coeftable[1,2])
  CtrySubsetResults[i,4] <- as.numeric(fereg$coeftable[1,3])
  CtrySubsetResults[i,5] <- as.numeric(fereg$coeftable[1,4])
}

CtrySubsetResults %>%
  ggplot() +
  geom_histogram(aes(x = Coef), binwidth = 0.0001) +
  geom_vline(xintercept = mean(CtrySubsetResults$Coef), linetype = 'dashed', alpha = 0.5) +
  theme_minimal() +
  coord_cartesian(xlim = c(0.005, 0.01))

#### Table A1 ####

OlsData %>%
  distinct(ccode) %>%
  mutate(`Country Name` = countrycode(as.numeric(ccode), 'cown', 'country.name')) %>%
  arrange(`Country Name`) 

#### Table A2 ####

OlsData %>% 
  dplyr::select(uri_p, sri_p, eu_stock, usa_stock, ofdi_stock_context_index, ofdi_stock_context_uri_p, ofdi_stock_context_sri_p, v2x_frassoc_thick, gsp, gdp_growth, UrbanPop, ilo_rat, gdppc_wanlin, export_devd_gdp, TotalOfdi, conflict, io_mem_n, ngoc) %>%
  mutate(across(uri_p:ngoc, \(x) as.numeric(x))) %>%
  summarise(across(.cols = c(uri_p, sri_p, eu_stock, usa_stock, ofdi_stock_context_index, ofdi_stock_context_uri_p, ofdi_stock_context_sri_p, v2x_frassoc_thick, gsp, gdp_growth, UrbanPop, ilo_rat, gdppc_wanlin, export_devd_gdp, TotalOfdi, conflict, io_mem_n, ngoc),
                   .fns = list('Min' = min,'Mean' = mean, 'Median' = median, 'Max' = max,'SD' = sd),
                   na.rm = T)) %>%
  pivot_longer('uri_p_Min':'ngoc_SD', names_to = 'Variable', values_to = ' ') %>%
  mutate(Stat = str_extract_all(Variable, 'Mean|Median|SD|Min|Max'),
         Variable = str_remove_all(Variable, '_Mean|_Median|_SD|_Min|_Max')) %>%
  pivot_wider(id_cols = 'Variable', names_from = 'Stat', values_from = ' ') %>%
  mutate(Source = c(rep('Barry et al. (2022)', 2),
                    rep('UNCTAD Bilateral FDI Database', 2),
                    'UNCTAD, OECD/AIAS ICTWSS & ILOSTAT',
                    rep('UNCTAD, Barry et al. (2022)',2),
                    'Coppedge et al. (2022)',
                    'Sari et al. (2016)',
                    'ILOSTAT',
                    rep('World Bank WDI', 3),
                    'IMF DOTS',
                    'UNCTAD Bilateral FDI Database',
                    'PRIO',
                    'Schofer et al. (2022)',
                    'Rocabert (2021)'),
         Variable = c('Union Rights (URI)',
                 'Substantial Rights (SRI)',
                 'Stock FDI in Europe',
                 'Stock FDI in USA',
                 'OFDI Context - Institutional',
                 'OFDI Context - Workr URI',
                 'OFDI Context - Workr SRI',
                 'Freedom of Political Association (V-Dem)',
                 'GSP',
                 'Total ILO Ratification',
                 'GDP Growth',
                 'Urban Population (%)',
                 'GDP per capita (constant 2010 USD)',
                 'Exports to Developed (% GDP)',
                 'Total OFDI Stock',
                 'Conflict - PRIO',
                 'INGO Presence',
                 'IO Membership'))

#### Table A3 ####

screenreg(Tables1andA3[c(9:12)],
        digits = 3,
        custom.header = list('FDI to EU' = 1:2, 'FDI to USA' = 3:4),
        custom.model.names = c('URI',
                               'SRI',
                               'URI',
                               'SRI'),
        custom.coef.map = list('eu_stock_lag1'='OFDI Stock t-1',
                               'usa_stock_lag1'='OFDI Stock t-1',
                               'tfp_lag1'='Total Factor Productivity t-1',
                               'plp_lag1'='Potential Labor Power t-1',
                               'indva_lag1'='Industry VA',
                               'gov1rlc_lag1'='Legislative Ideology',
                               'log(1 + InwardFDIstock_lag1)'='Inward FDI Stock t-1',
                               'v2x_frassoc_thick_lag1'='Freedom of Assoc. t-1',
                               'gsp_lag1'='GSP t-1',
                               'log(1 + ilo_rat_lag1)'='ILO Ratification t-1',
                               'log(1 + gdppc_wanlin_lag1)'='GDP per capita t-1',
                               'gdp_growth_lag1'= 'GDP Growth t-1',
                               'UrbanPop_lag1'='Urban Population t-1',
                               'log(1 + export_devd_gdp_lag1)'='Exports to Devd t-1',
                               'log(1 + TotalOfdi_lag1)' = 'Total OFDI Stock',
                               'conflict1' = 'Conflict',
                               'log(1 + io_mem_n)' = 'IO Membership',
                               'log(1 + ngoc)' = 'INGO Presence'),
        custom.gof.rows = list('Country Fixed Effects' = c('Yes', 'Yes', 'Yes', 'Yes'),
                               'Year Fixed Effects' = c('Yes', 'Yes', 'Yes', 'Yes')),
        stars = c(0.01, 0.05, 0.1))

#### Table A4 ####

CtrySubsetResults %<>%
  arrange(OmittedCountry)

Encoding(CtrySubsetResults$OmittedCountry) <- 'latin1'

data.table(format(CtrySubsetResults, digits = 3))

#### Table A5 ####

ivs_vec <- c('eu_stock_lag1',
             'usa_stock_lag1')
dvs_vec <- c('KS_facb',
             'marx_facb')

ivs <- paste0(" ~ ", ivs_vec)

dvs_ivs <- c(unlist(lapply(ivs, function(x) paste0(dvs_vec, x,
                                                   '+
                v2x_frassoc_thick_lag1 +
               gsp_lag1 + 
               log(1 + ilo_rat_lag1) +
               log(1 + gdppc_wanlin_lag1)+
               gdp_growth_lag1 +
               UrbanPop_lag1+
               log(1 + export_devd_gdp_lag1) + 
               conflict +
               log(1 + TotalOfdi_lag1) +
               log(1 + io_mem_n) + 
               log(1 + ngoc)|ccode+year'))))
formulas <- lapply(dvs_ivs, formula)

TableA5 <- lapply(formulas, function(x) {
  feols(x, data = OlsData,
        se = 'cluster')
})

screenreg(TableA5,
          digits = 3,
          custom.header = list('FDI to EU' = 1:2, 'FDI to USA' = 3:4),
          custom.model.names = c('Kucer-Sari',
                                 'Marx et al.',
                                 'Kucer-Sari',
                                 'Marx et al.'),
          custom.coef.map = list('eu_stock_lag1'='OFDI Stock t-1',
                                 'usa_stock_lag1'='OFDI Stock t-1'),
          custom.gof.rows = list('Country Fixed Effects' = c('Yes', 'Yes', 'Yes', 'Yes'),
                                 'Year Fixed Effects' = c('Yes', 'Yes', 'Yes', 'Yes')),
          stars = c(0.01, 0.05, 0.1)
)

#### Table A6 ####

ivs_vec <- c('eu_stock_lag1',
             'usa_stock_lag1')
dvs_vec <- c('uri_l',
             'sri_l')

ivs <- paste0(" ~ ", ivs_vec)

dvs_ivs <- c(unlist(lapply(ivs, function(x) paste0(dvs_vec, x,
                                                   '+
               v2x_frassoc_thick_lag1 +
               gsp_lag1 + 
               gdp_growth_lag1 +
               UrbanPop_lag1+
               log(1 + ilo_rat_lag1) +
               log(1 + gdppc_wanlin_lag1)+
               log(1 + export_devd_gdp_lag1) + 
               log(1 + InwardFDIstock_lag1) +
               conflict +
               log(1 + TotalOfdi_lag1) +
               log(1 + io_mem_n) + 
               log(1 + ngoc)|ccode+year'))))
formulas <- lapply(dvs_ivs, formula)

TableA6 <- lapply(formulas, function(x) {
  feols(x, data = OlsData,
        se = 'cluster')
})

screenreg(TableA6,
        digits = 3,
        custom.header = list('FDI to EU' = 1:2, 'FDI to USA' = 3:4),
        custom.model.names = rep(c('URI',
                                   'SRI'),2),
        custom.coef.map = list('eu_stock_lag1'='OFDI Stock t-1',
                               'usa_stock_lag1'='OFDI Stock t-1'),
        custom.gof.rows = list('All Controls' = c('Yes', 'Yes', 'Yes', 'Yes'),
                               'Country Fixed Effects' = c('Yes','Yes', 'Yes', 'Yes'),
                               'Year Fixed Effects' = c( 'Yes', 'Yes', 'Yes', 'Yes')),
        stars = c(0.01, 0.05, 0.1))


